library(sf)
library(readxl)
library(dplyr)
library(plyr)
library(ggplot2)
library(afrihealthsites)
library(ggpubr)
library(afriadmin)
library(tmap)
library(cowplot)

# Install Malawi MFL
malawi_MFL = read_excel("~/malawi-health-facilities-1/MHFR_Facilities 1.xlsx")

# Convert to sf

## omit NA's
new_malawi_MFL = na.omit(malawi_MFL)

## check for NA 
any(is.na(new_malawi_MFL))
## [1] FALSE
## transform geometry columns into numeric 
sapply(new_malawi_MFL, class)
##        CODE        NAME COMMON NAME   OWNERSHIP        TYPE      STATUS 
## "character" "character" "character" "character" "character" "character" 
##        ZONE    DISTRICT DATE OPENED    LATITUDE   LONGITUDE 
## "character" "character" "character" "character" "character"
new_malawi_MFL = transform(new_malawi_MFL, LATITUDE = as.numeric(LATITUDE), 
                                           LONGITUDE = as.numeric(LONGITUDE))
## Warning in eval(substitute(list(...)), `_data`, parent.frame()): NAs introduced
## by coercion
any(is.na(new_malawi_MFL)) ## check for NA
## [1] TRUE
new_malawi_MFL = na.omit(new_malawi_MFL) ## and omit

## convert to sf object
malawi_facilities_MFL = st_as_sf(new_malawi_MFL, coords = c("LONGITUDE", "LATITUDE"), dim = "XY")

malawi_facilities_MFL = st_set_crs(malawi_facilities_MFL, 4326) ## set CRS, is WGS84 right?

head(malawi_facilities_MFL)
## Simple feature collection with 6 features and 9 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 33.74129 ymin: -15.84 xmax: 35.09 ymax: -13.79742
## geographic CRS: WGS 84
##       CODE                                               NAME  COMMON.NAME
## 1 MC010002                               A + A private clinic          A+A
## 2 BT240003                                       A-C Opticals A.C Opticals
## 3 BT240005                                Akwezeke PVT Clinic Akwezeke Pvt
## 4 BT240006                                  AB Medical Clinic        Abowa
## 5 LL040007                                 ABC Comm. Hospital   ABC Clinic
## 6 LL040010 Achikondi Women Community Friendly Services Clinic    Achikondi
##                                       OWNERSHIP       TYPE     STATUS
## 1                                       Private     Clinic Functional
## 2                                       Private     Clinic Functional
## 3                                       Private     Clinic Functional
## 4                                       Private     Clinic Functional
## 5 Christian Health Association of Malawi (CHAM)   Hospital Functional
## 6                                       Private Dispensary Functional
##                 ZONE DISTRICT DATE.OPENED                   geometry
## 1 Centrals West Zone  Mchinji  Jan 1st 75 POINT (33.88563 -13.79742)
## 2    South East Zone Blantyre  Jan 1st 75        POINT (35.03 -15.8)
## 3    South East Zone Blantyre  Jan 1st 75       POINT (35.09 -15.84)
## 4    South East Zone Blantyre  Jan 1st 75       POINT (35.09 -15.84)
## 5 Centrals West Zone Lilongwe  Jan 1st 75 POINT (33.74129 -13.96816)
## 6 Centrals West Zone Lilongwe  Jan 1st 75  POINT (33.7793 -13.95473)

Overview:

  1. Malawi MFL
# Re-order the facility types
malawi_MFL$TYPE = as.factor(malawi_MFL$TYPE)

malawi_MFL$TYPE = factor(malawi_MFL$TYPE, levels = c("Central Hospital", "District Hospital", "Hospital", "Health Centre", "Clinic", "Health Post", "Dispensary", "Private", "Unclassified"))

# Number of each type + ownership
facility_types_MFL = as.data.frame(table(malawi_MFL$TYPE, malawi_MFL$OWNERSHIP))
head(facility_types_MFL)
##                Var1            Var2 Freq
## 1  Central Hospital Aquaid Lifeline    0
## 2 District Hospital Aquaid Lifeline    0
## 3          Hospital Aquaid Lifeline    0
## 4     Health Centre Aquaid Lifeline    0
## 5            Clinic Aquaid Lifeline    0
## 6       Health Post Aquaid Lifeline    0
## bar plot of no. of facility types 
plot_facility_types_MFL = ggplot(facility_types_MFL, aes(x=Var1, y=Freq, fill=Var2)) + geom_bar(position = "stack", stat = "identity")
plot_facility_types_MFL = plot_facility_types_MFL + labs(x = "Facility types", y = "Frequency", fill="Ownership") + scale_fill_brewer(palette = "Set2") + coord_flip() + theme_pubclean() + ggtitle("MFL") + theme(legend.title = element_text(size = 9), legend.text = element_text(size = 9), plot.title = element_text(face = "bold", hjust = 0), legend.position = "bottom")
plot_facility_types_MFL

par(mar=c(11,4,4,4))
# Re-order ownership
malawi_MFL$OWNERSHIP = as.factor(malawi_MFL$OWNERSHIP)

malawi_MFL$OWNERSHIP = factor(malawi_MFL$OWNERSHIP, levels = c("Government", "Private", "Christian Health Association of Malawi (CHAM)", "Non-Government", "Mission/Faith-based (other than CHAM)", "Other", "Parastatal", "Aquaid Lifeline"))

# Number of each type of ownership
ownership_MFL = as.data.frame(table(malawi_MFL$OWNERSHIP))
ownership_MFL
##                                            Var1 Freq
## 1                                    Government  695
## 2                                       Private  495
## 3 Christian Health Association of Malawi (CHAM)  192
## 4                                Non-Government   69
## 5         Mission/Faith-based (other than CHAM)   62
## 6                                         Other   27
## 7                                    Parastatal    5
## 8                               Aquaid Lifeline    1
## bar plot of ownership
plot_ownership_MFL = ggplot(ownership_MFL, aes(x=Var1, y=Freq)) + geom_bar(stat = "identity", fill="slategray")
plot_ownership_MFL = plot_ownership_MFL + labs(y = "Frequency") + coord_flip() + theme_pubclean() + ggtitle("MFL") + theme(axis.title.y = element_blank(), plot.title = element_text(face = "bold"))
plot_ownership_MFL

  1. WHO-KWTRP
# Malawi WHO data.frame
malawi_WHO <- afrihealthsites("malawi", datasource='who', plot=FALSE, returnclass='dataframe')

head(malawi_WHO)
## # A tibble: 6 x 10
##   Country Admin1 `Facility name` `Facility type` Ownership   Lat  Long
##   <chr>   <chr>  <chr>           <chr>           <chr>     <dbl> <dbl>
## 1 Malawi  Centr… 80 Block Clinic Clinic          MoH       -12.9  33.4
## 2 Malawi  Centr… ABC Community … Clinic          FBO       -14.0  33.7
## 3 Malawi  Centr… Adventist Heal… Health Centre   FBO       -14.0  33.8
## 4 Malawi  Centr… Alinafe Commun… Community Hosp… FBO       -13.4  34.2
## 5 Malawi  Centr… Area 18 Health… Health Centre   MoH       -13.9  33.8
## 6 Malawi  Centr… Area 25 Health… Health Centre   MoH       -13.9  33.8
## # … with 3 more variables: `LL source` <chr>, iso3c <chr>,
## #   facility_type_9 <chr>
# Re order facility types and ownership 
malawi_WHO$`Facility type` = as.factor(malawi_WHO$`Facility type`)
malawi_WHO$`Facility type` = factor(malawi_WHO$`Facility type`, levels = c("Central Hospital", "District Hospital", "Mission Hospital", "Rural Hospital", "Community Hospital", "Health Centre", "Clinic", "Health Post/Dispensary"))

malawi_WHO$Ownership = as.factor(malawi_WHO$Ownership)
malawi_WHO$Ownership = factor(malawi_WHO$Ownership, levels = c("MoH", "FBO", "Local authority", "NGO"))

# No. of original facility types + ownership
facility_types_WHO = as.data.frame(table(malawi_WHO$`Facility type`, malawi_WHO$Ownership))
head(facility_types_WHO)
##                 Var1 Var2 Freq
## 1   Central Hospital  MoH    4
## 2  District Hospital  MoH   24
## 3   Mission Hospital  MoH    0
## 4     Rural Hospital  MoH   17
## 5 Community Hospital  MoH    0
## 6      Health Centre  MoH  332
## bar plot of original facility types 
plot_facility_types_WHO = ggplot(facility_types_WHO, aes(x=Var1, y=Freq, fill=Var2)) + geom_bar(position = "stack", stat = "identity")
plot_facility_types_WHO = plot_facility_types_WHO + labs(x = "Facility types", y = "Frequency", fill="Ownership") + scale_fill_brewer(palette = "Set2") + coord_flip() + theme_pubclean() + ggtitle("WHO") + theme(legend.title = element_text(size = 9), legend.text = element_text(size = 9), plot.title = element_text(face = "bold", hjust = 0), legend.position = "bottom") + expand_limits(y=c(0,600))
plot_facility_types_WHO

# Re order
malawi_WHO$facility_type_9 = as.factor(malawi_WHO$facility_type_9)
malawi_WHO$facility_type_9 = factor(malawi_WHO$facility_type_9, levels = c("Hospital", "Health Centre", "Health Clinic", "Health Post", "Community Health Unit"))

# No. of reclassified facility types
RC_facility_types_WHO = as.data.frame(table(malawi_WHO$facility_type_9))
RC_facility_types_WHO
##                    Var1 Freq
## 1              Hospital   80
## 2         Health Centre  457
## 3         Health Clinic   22
## 4           Health Post   87
## 5 Community Health Unit    2
## bar plot of reclassified facility types 
plot_RC_facility_types_WHO = ggplot(RC_facility_types_WHO, aes(x=Var1, y=Freq)) + geom_bar(stat = "identity", fill="slategray")
plot_RC_facility_types_WHO = plot_RC_facility_types_WHO + labs(y = "Frequency") + coord_flip() + theme_pubclean() + theme(axis.title.y = element_blank()) + ggtitle("Reclassified facility types")
plot_RC_facility_types_WHO

# Types of ownership
ownership_WHO = as.data.frame(table(malawi_WHO$Ownership))
ownership_WHO
##              Var1 Freq
## 1             MoH  467
## 2             FBO  173
## 3 Local authority    5
## 4             NGO    3
## bar plot of ownership
plot_ownership_WHO = ggplot(ownership_WHO, aes(x=Var1, y=Freq)) + geom_bar(stat = "identity", fill="slategray")
plot_ownership_WHO = plot_ownership_WHO + labs(y = "Frequency") + coord_flip() + theme_pubclean() + ggtitle("WHO") + theme(axis.title.y = element_blank(), plot.title = element_text(face = "bold")) + expand_limits(y=c(0,600))
plot_ownership_WHO

Both data sources contain no information on services available, capacity or equipment. MFL does state whether facility is functional.

Classification of MFL facilities aligns more with the structure of the health care system in Malawi (community, primary, secondary, tertiary), it differentiates central hospitals from district and other hospitals. WHO has additional rural and mission hospitals, where do they fit in?

https://www.health.gov.mw/index.php/2016-01-06-19-58-23/national-aids states that at community level, health posts, dispensaries and maternity clinics offer services. Primary includes health centers and community hospitals, secondary consists of district and some CHAM hospitals, tertiary includes central hospitals.

Analysis:

## tmap mode set to interactive viewing
## Warning in sf::st_is_longlat(shp2): bounding box has potentially an invalid
## value range for longlat data
## Warning in sf::st_is_longlat(shp2): bounding box has potentially an invalid
## value range for longlat data
## Warning in sf::st_is_longlat(shp2): bounding box has potentially an invalid
## value range for longlat data

Qs to address?:

  1. How many facilities are in the same location across the MFL and WHO datasets?
  2. Do same facilities share same names/other attributes?
  3. The ones that aren’t, are they within 50m of another facility?
# Qs 1 - how many intersect?

## convert malawi_WHO to sf object

class(malawi_WHO)
## [1] "tbl_df"     "tbl"        "data.frame"
any(is.na(malawi_WHO))
## [1] TRUE
new_malawi_WHO = na.omit(malawi_WHO) ## omit NA

sf_malawi_WHO = st_as_sf(new_malawi_WHO, coords = c("Long", "Lat"), dim = "XY")
sf_malawi_WHO = st_set_crs(sf_malawi_WHO, 4326)

## st_intersection
intersect_WHO_MFL = st_intersection(x=sf_malawi_WHO, y=malawi_facilities_MFL)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
intersect_WHO_MFL ## only 2 intersect directly, so are same up to 5 decimal places?
## Simple feature collection with 2 features and 17 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 33.29456 ymin: -11.53894 xmax: 33.41925 ymax: -11.45836
## geographic CRS: WGS 84
## # A tibble: 2 x 18
##   Country Admin1 Facility.name Facility.type Ownership LL.source iso3c
## * <chr>   <chr>  <chr>         <fct>         <fct>     <chr>     <chr>
## 1 Malawi  North… Euthini Heal… Health Centre MoH       GPS       MWI  
## 2 Malawi  North… Madede Healt… Health Centre MoH       GPS       MWI  
## # … with 11 more variables: facility_type_9 <fct>, CODE <chr>, NAME <chr>,
## #   COMMON.NAME <chr>, OWNERSHIP <chr>, TYPE <fct>, STATUS <chr>, ZONE <chr>,
## #   DISTRICT <chr>, DATE.OPENED <chr>, geometry <POINT [°]>
# Qs 2 - Do they share same attributes
## ownership is same, 1 name is same, Euthini registered as a hospital in MFL but as a health centre in WHO

# Qs 3 - how many/which are within 50m of another facility?